#!/usr/bin/perl -w
use strict;

use FindBin;
use lib ("$FindBin::Bin/..", "/net/cpp-group/Leo/bin");
use db_parameters;
use ortho_species;
use Spreadsheet::WriteExcel::Big;
use Data::Dump qw(dump);
use dbh_helper;
use DBI;
use xls_formats;

print STDERR <<"HEADLINE";
888888888888888888888888888888888888888888888888888888888888888888888

    xls_protocols_counts

		gene counts for each stage in each protocol

888888888888888888888888888888888888888888888888888888888888888888888

HEADLINE

#_________________________________________________________________________________________
#
#	get gene counts
#
#_________________________________________________________________________________________
sub get_gene_counts_by_category($\%)
{
	print STDERR "\tGet gene counts by category from panda...\n";
	my ($dbh, $gene_counts) = @_;
    my %counts;
	my $sql_cmd = <<"PL/SQLCMD";
	SELECT 
		protocols.name,
		species,
		gene_type, 
		count
	FROM
		orthologs.predicted_type_summary NATURAL JOIN 
		orthologs.protocols 
	ORDER BY 
		protocols.name, 
		species, 
		gene_type;
PL/SQLCMD
	print STDERR "\tGet counts from panda...\n";
	my $data = $dbh->selectall_arrayref($sql_cmd) or die $dbh->errstr.".";
	for my $line(@$data)
	{
		my ($protocol_name, $species, $gene_type, $cnt) = @$line;
		$gene_counts->{$protocol_name}{$species}{$gene_type} = $cnt;
	}

}



#_________________________________________________________________________________________
#
#	get dupications
#
#_________________________________________________________________________________________
sub get_duplications_count($\%)
{
	print STDERR "\tGet duplications from panda...\n";
	my ($dbh, $cnt_genes) = @_;
	my $sql_cmd = <<"PL/SQLCMD";
    SELECT
        protocol_name,
        species,
        SUM(cnt_duplicates) 
    FROM 
    (
        SELECT 
            COUNT(*) - 1 AS cnt_duplicates,
            ortholog_id,
            species,
            protocols.name AS protocol_name
            
        FROM
            orthologs.ortholog_sets NATURAL JOIN 
            orthologs.protocols
        WHERE
            (ortholog_type ~ 'to' OR
             ortholog_type ~ 'inparalog')
        GROUP BY 
            ortholog_id,
            protocols.name,
            species
    ) as a
    GROUP BY 
        protocol_name,
        species;
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die $dbh->errstr.".";
	for my $line(@$data)
	{
		my ($protocol_name, $species, $cnt) = @$line;
		$cnt_genes->{DUPLICATIONS}{$protocol_name}{$species} = $cnt;
	}

}


#_________________________________________________________________________________________
#
#	get original gene count
#
#_________________________________________________________________________________________
sub get_original_gene_count($\%)
{
	print STDERR "\tGet original gene counts from panda...\n";
	my ($dbh, $cnt_genes) = @_;
	my $sql_cmd = <<"PL/SQLCMD";
    SELECT 
            protocols.name, 
            predicted_type.species, 
            COUNT(DISTINCT gene_id),
            COUNT(DISTINCT prot_id) AS COUNT 
        FROM 
            orthologs.predicted_type JOIN 
            ens_id USING (gene_id) NATURAL JOIN 
            orthologs.protocols 
        GROUP BY 
            protocols.name, 
            predicted_type.species 
        ORDER BY 
            name, 
            predicted_type.species;
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die $dbh->errstr.".";
	for my $line(@$data)
	{
		my ($protocol_name, $species, $cnt_g, $cnt_p) = @$line;
		$cnt_genes->{ORIGINAL_GENES}{$protocol_name}{$species} = $cnt_g;
		$cnt_genes->{ORIGINAL_TRANSCRIPTS}{$protocol_name}{$species} = $cnt_p;
	}
}


#_________________________________________________________________________________________
#
#	get original gene count
#
#_________________________________________________________________________________________
sub get_split_genes_count($\%)
{
	print STDERR "\tGet split/demerged gene counts from panda...\n";
	my ($dbh, $cnt_genes) = @_;
	my $sql_cmd = <<"PL/SQLCMD";
    SELECT 
            protocols.name, 
            species, 
            COUNT(*), 
            COUNT(DISTINCT gene_id) 
        FROM
            split_genes NATURAL JOIN 
            ens_id NATURAL JOIN 
            protocols 
        GROUP BY 
            species, name; 
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die $dbh->errstr.".";
	for my $line(@$data)
	{
		my ($protocol_name, $species, $cnt_i, $cnt_f) = @$line;
		$cnt_genes->{SPLIT_INTO}{$protocol_name}{$species} = $cnt_i;
		$cnt_genes->{SPLIT_FROM}{$protocol_name}{$species} = $cnt_f;
	}
}



#-----------------------------------------------------------------------------------------

#	print_gene_counts_diagram

#-----------------------------------------------------------------------------------------
sub print_gene_counts_diagram($$\%\%\%)
{
	my ($dbh, $workbook, $formats, $counts,
		$extra_gene_counts) = @_;
	my @protocol_names = ("homo_mono_ens", "Andreas Monodelphis 2 Gene Predictions", "Andreas Monodelphis Paralogue Repredictions");
	my $worksheet = $workbook->addworksheet("Gene counts");
	my @species = ($ortho_name1, $ortho_name2);
	my %column_names = (
					"1 to 1"									=>  "1 to 1", 
					"m to 1"									=>  "m to 1", 
					"1 to m"									=>  "1 to m", 
					"m to m"									=>  "m to m", 
					"Orphan inparalog inconsistent phylogeny"	=>  "Duplicated with inconsistent phylogeny", 
					"Orphan singleton inconsistent phylogeny"	=>  "Singletons with inconsistent phylogeny", 
					"Orphan inparalog"							=>  "Orphaned inparalogues", 
					"Orphan low dS homolog other species"		=>  "Orphans with close [OTHER] homologues", 
					"Orphan low dS homolog within species"		=>  "Orphans with close [SAME] homologues", 
					"Orphan homolog other species"				=>  "Orphans with only distant [OTHER] homologues", 
					"Orphan homolog within species"				=>  "Orphans with only distant [SAME] homologues", 
					"Orphan No homolog"							=>  "No homology (e < 10-5, 50 codons)",
					"pseudo"									=>  "Pseudogenes", 
					"Orphan by pseudo"							=>  "Orthologues of pseudogenes");
	my @categories = (  "1 to 1", "m to 1", "1 to m", "m to m",
					"Orphan inparalog inconsistent phylogeny",
					"Orphan singleton inconsistent phylogeny",
					"Orphan inparalog",
					"Orphan low dS homolog other species", 
					"Orphan low dS homolog within species", 
					"Orphan homolog other species",
					"Orphan homolog within species",
					"Orphan No homolog", 
					"pseudo", "Orphan by pseudo");


	print STDERR "\tPrint Orthology counts to Excel...\n";


	#
	#	write column names
	#
	{
		$worksheet->set_column(1, 1, 50);
		# extra line for species name
		my $row = 3;
		for my $species_id(0,1)
		{
			my $species			= $species[$species_id];
			my $other_species	= $species[($species_id + 1) % 2];

			#
			# PhyOP subtype
			#
			for my $item_name(@categories)
			{
				my $descriptive_name = $column_names{$item_name};
				$descriptive_name =~ s/\[SAME\]/$species/;
				$descriptive_name =~ s/\[OTHER\]/$other_species/;
				$worksheet->write($row++,1, $descriptive_name, $formats->{blue_bold});
			}

			#spacer
			++$row;

			#
			# Other gene count categories
			#
			for my $category ("All PhyOP Genes", "Original Genes", "Original Transcripts",
							  "Duplications", "Demerged Genes", "Split from")
			{
				$worksheet->write($row++,1, $category, $formats->{blue_bold});
			}



			# spacer
			# extra line for species name
			$row += 3;
		}
	}
	
	#
	#	write data
	#
	my $col = 2;
	for my $protocol(@protocol_names)
	{
		# set column width and write protocol_name
		$worksheet->set_column($col, $col, 20);
		$worksheet->write(1, $col, $protocol, $formats->{red_bold_yellow});
		
		my $row = 2;
		for my $species(@species)
		{
			$worksheet->write($row++, $col, $species, $formats->{blue_bold});
			unless (exists $counts->{$protocol})
			{
				warn "Error:\tData for $protocol does not exist\n";
			}
			
			for my $item_name(@categories)
			{
				$extra_gene_counts->{TOTAL}{$protocol}{$species} 
					+= $counts->{$protocol}{$species}{$item_name}
						if (exists $counts->{$protocol}{$species}{$item_name});
			}

			for my $item_name(@categories)
			{
				$worksheet->write($row, $col, 
								  $counts->{$protocol}{$species}{$item_name}, 
								  $formats->{comma})
					if (exists $counts->{$protocol}{$species}{$item_name});
				$row++;
			}
			

			#spacer
			++$row;

			#original counts, duplications, and merged gene counts
			for my $category (qw(TOTAL ORIGINAL_GENES ORIGINAL_TRANSCRIPTS DUPLICATIONS 
								 SPLIT_INTO SPLIT_FROM) )
			{
				$worksheet->write($row, $col, 
								  $extra_gene_counts->{$category}{$protocol}{$species}, 
								  $formats->{comma})
					if (exists $extra_gene_counts->{$category}{$protocol}{$species});
				++$row;

			}

		
			$row+=2;
			
		}
		$col++;
	}
}


#-----------------------------------------------------------------------------------------
#	Get DB parameters
#
use db_parameters;
#my @db_param = get_db_parameters();
#
#	Get DB parameters
#-----------------------------------------------------------------------------------------






#
# connect to genome database
#
my $dbh = connect_to_panda();

	#
	# get counts by category
	#
    my %counts;
	get_gene_counts_by_category($dbh, %counts);


	#
	# Get other counts
	#
	my %extra_counts;

	# get original gene counts
	get_original_gene_count($dbh, %extra_counts);

	# get original gene counts
	get_duplications_count($dbh, %extra_counts);

	# get split/demerged gene counts
	get_split_genes_count($dbh, %extra_counts);



	
	#
	#	workbook
	#
	my $workbook  = Spreadsheet::WriteExcel::Big->new("$dir_orthologs_output/excel/protocols_gene_counts.xls");
	die "Problems creating new Excel file: $!" unless defined $workbook;

	my %formats;
	get_xls_formats($workbook, %formats);


	print_gene_counts_diagram($dbh, $workbook, %formats, %counts, 
							  %extra_counts);
	$workbook->close();




print STDERR "\tCompleted\n";



$dbh->disconnect();





